function [delays,param] = txdelay(varargin)

%TXDELAY   Transmit delays for a linear or convex array
%   TXDELAY returns the transmit time delays for focused, plane or circular
%   beam patterns in a linear or convex array.
%
%   DELAYS = TXDELAY(X0,Z0,PARAM) returns the transmit time delays which
%   must be used to generate a pressure field focused at the point
%   (x0,z0). Note: If z0 is negative, then the point (x0,z0) is a virtual
%   source. The properties of the medium and linear array must be given in
%   the structure PARAM (see below).
%
%   DELAYS = TXDELAY(PARAM,TILT) returns the transmit time delays which
%   must be used to get a tilted plane wave. TILT is the tilt angle, i.e
%   the angle between the emission axis and Z-vertical axis.
%
%   DELAYS = TXDELAY(PARAM,TILT,WIDTH) yields the transmit time delays
%   necessary for creating a circular wave. The sector enclosed by the
%   circular waves is characterized by the angular width and tilt. TILT
%   represents the sector tilt, WIDTH is the sector width (both in
%   radians). This option is not available for a convex array.
%
%   X0, Z0, TILT and WIDTH can be vectors. In that case, DELAYS is a matrix
%   whose rows contain the different delay laws.
%
%   [DELAYS,PARAM] =  TXDELAY(...) updates the PARAM structure parameters
%   including the default values. PARAM will also include PARAM.TXdelay
%   which is equal to DELAYS (in s).
%
%   [...] = TXDELAY (no input parameter) runs an interactive example
%   simulating a focused pressure field generated by a 2.7 MHz phased
%   array. The user must choose the focus position.
%
%   Units: X0,Z0 must be in m; TILT, WIDTH must be in rad. DELAYS are in s.
%
%   PARAM is a structure which must contain the following fields:
%   ------------------------------------------------------------
%   1) PARAM.pitch: pitch of the linear array (in m, REQUIRED)
%   2) PARAM.Nelements: number of elements in the transducer array (REQUIRED)
%   3) PARAM.radius: radius of curvature (in m, default = Inf)
%   4) PARAM.c: longitudinal velocity (in m/s, default = 1540 m/s)
%
%   ---
%   NOTE #1: X- and Z-axes
%   The X axis is PARALLEL to the transducer and points from the first
%   (leftmost) element to the last (rightmost) element (X = 0 at the CENTER
%   of the transducer). The Z axis is PERPENDICULAR to the transducer and
%   points downward (Z = 0 at the level of the transducer, Z increases as
%   depth increases).
%   ---
%   NOTE #2: TILT describes the tilt angle in the trigonometric direction.
%   ---
%
%   Example #1:
%   ----------
%   %-- Generate a focused pressure field with a phased-array transducer
%   % Phased-array @ 2.7 MHz:
%   param = getparam('P4-2v');
%   % Focus position:
%   x0 = 2e-2; z0 = 5e-2;
%   % TX time delays:
%   dels = txdelay(x0,z0,param);
%   % Grid:
%   x = linspace(-4e-2,4e-2,200);
%   z = linspace(0,10e-2,200);
%   [x,z] = meshgrid(x,z);
%   % RMS pressure field:
%   P = pfield(x,z,dels,param);
%   imagesc(x(1,:)*1e2,z(:,1)*1e2,20*log10(P/max(P(:))))
%   hold on, plot(x0*1e2,z0*1e2,'k*'), hold off
%   colormap hot, axis equal tight ij
%   caxis([-20 0])
%   c = colorbar;
%   c.YTickLabel{end} = '0 dB';
%   xlabel('[cm]')
%
%   Example #2:
%   ----------
%   %-- Generate a plane wave a convex transducer
%   % Convex array @ 3.6 MHz:
%   param = getparam('C5-2v');
%   % Tilt angle = 10 degrees:
%   tilt = pi/18; % in rad
%   % TX apodization
%   param.TXapodization = [ones(1,100) zeros(1,28)];
%   % TX time delays:
%   dels = txdelay(param,tilt);
%   % 8cm-by-8cm grid:
%   [x,z] = impolgrid(100,8e-2,param);
%   % RMS pressure field:
%   P = pfield(x,z,dels,param);
%   pcolor(x*1e2,z*1e2,20*log10(P/max(P(:))))
%   shading interp
%   colormap hot, axis equal tight ij
%   caxis([-20 0])
%   c = colorbar;
%   c.YTickLabel{end} = '0 dB';
%   xlabel('[cm]')
%
%
%   This function is part of <a
%   href="matlab:web('https://www.biomecardio.com/MUST')">MUST</a> (Matlab UltraSound Toolbox).
%   MUST (c) 2020 Damien Garcia, LGPL-3.0-or-later
%
%   See also PFIELD, SIMUS, DAS, DASMTX, GETPARAM.
%
%   -- Damien Garcia -- 2015/03, last update: 2020/06/11
%   website: <a
%   href="matlab:web('https://www.biomecardio.com')">www.BiomeCardio.com</a>

%-- Check the input arguments
if nargin==0
    if nargout>0
        [delays,param] = RunTheExample;
    else
        RunTheExample;
    end
    return
elseif nargin==2 % Plane wave: TXDELAY(param,tilt)
    param = varargin{1};
    option = 'Plane Wave';
elseif nargin==3
    if isstruct(varargin{3}) % Origo: TXDELAY(x0,z0,param)
        param = varargin{3};
        option = 'Origo';
    else  % Circular wave: TXDELAY(param,tilt,width)
        param = varargin{1};
        option = 'Circular Wave';
    end
else
    error('Wrong input arguments.')
end
assert(isstruct(param),'Wrong input arguments. PARAM must be a structure.')

%-- Number of elements
if isfield(param,'Nelements')
    N = param.Nelements;
else
    error('The number of elements (PARAM.Nelements) is required.')
end

%-- Pitch (in m)
if ~isfield(param,'pitch')
    error('A pitch value (PARAM.pitch) is required.')
end

%-- Longitudinal velocity (in m/s)
if ~isfield(param,'c')
    param.c = 1540;
end
c = param.c;

%-- Radius of curvature (in m)
% for a convex array
if ~isfield(param,'radius')
    param.radius = Inf; % default = linear array
end
R = param.radius;
isLINEAR = isinf(R);



%-- Positions of the transducer elements
if isLINEAR
    % we have a LINEAR ARRAY
    x = ((0:N-1)-(N-1)/2)*param.pitch;
    z = zeros(1,N);
else
    % we have a CONVEX ARRAY
    chord = 2*R*...
        sin(asin(param.pitch/2/R)*(N-1));
    h = sqrt(R^2-chord^2/4);
    THe = linspace(atan2(-chord/2,h),atan2(chord/2,h),N);
    z = R*cos(THe);
    x = R*sin(THe);
    % z = z-R;
    z = z-h;
end


switch option
    %-----
    case 'Plane Wave'
        tilt = varargin{2}(:);
        assert(all(abs(tilt)<pi/2),...
            'The tilt angles must verify |tilt| < pi/2')
        if isLINEAR
            delays = x.*sin(-tilt)/c;
        else
            % we have a CONVEX ARRAY
            
            % intersection point between the wavefront and the transducer
            xn = -R*sin(tilt); zn = R*cos(tilt)-h;

            % Note:
            % Equation of the line tangent to the transducer at (xn,zn):
            % X = -xn/(zn+h)*(X-xn) + zn
            
            % distances between this line and the elements (x,z)
            d = abs(z+xn./(zn+h)*x-xn.^2./(zn+h)-zn)./...
                sqrt(1+xn.^2./(zn+h).^2);
            
            delays = -d/c;
        end
    %-----
    case 'Origo'
        x0 = varargin{1}(:);
        z0 = varargin{2}(:);
        assert(numel(x0)==numel(z0),...
            'X0 and Z0 must have the same length.')
        delays = sqrt((x-x0).^2 + (z-z0).^2)/c;
        if isLINEAR
            delays = -delays.*sign(z0);
        elseif sqrt(x0.^2 + (R-z0).^2)<R
            delays = -delays;
        end
    %-----
    case 'Circular Wave'
        assert(isLINEAR,['Only the ''TXDELAY(X0,Z0,PARAM)'' syntax',...
            ' is available with a convex array.'])
        tilt = varargin{2}(:);
        width = varargin{3}(:);
        assert(numel(tilt)==numel(width),...
            'TILT and WIDTH must have the same length.')
        assert(all(width>0 & width<pi),...
            'The width angles must verify width > 0 and width < pi')
        L = (N-1)*param.pitch;
        %-- Origo
        [x0,z0] = angles2origo(L,tilt,width);
        %--
        delays = sqrt((x-x0).^2 + z0.^2)/c;
        delays = -delays.*sign(z0);
end

delays = delays-min(delays,[],2);

if nargout>1, param.TXdelay = delays; end

end

function [x0,z0] = angles2origo(L,tilt,width)
% Origo (virtual source) from the tilt and width angles
tilt = mod(tilt+pi/2,2*pi)-pi/2;
SignCorrection = ones(size(tilt));
idx = abs(tilt)>pi/2;
tilt(idx) = pi-tilt(idx);
SignCorrection(idx) = -1;
z0 = SignCorrection.*L./(tan(tilt-width/2)-tan(tilt+width/2));
x0 = SignCorrection.*z0.*tan(width/2-tilt)+L/2;
end

function [dels,param] = RunTheExample

% 2.7 MHz phased-array
param = getparam('P4-2v');
% Location of the focus
h = figure('Name','MUST - Matlab UltraSound Toolbox - Damien Garcia',...
    'NumberTitle','off');
plot((param.Nelements-1)*param.pitch/2*[-1 1],[0 0],'k','Linewidth',10)
axis equal ij
axis([-4e-2 4e-2 0 10e-2])
xlabel('x-position (m)')
ylabel('z-position (m)')
hp = helpdlg({'We have a 2.7 MHz phased array.',...
    'Choose the focus point in the figure.'},'Focus point Selection');
waitfor(hp)
[x0,z0] = ginput(1);
hold on
plot(x0,z0,'k*')
pause(0.2)
close(h)
% TX time delays:
dels = txdelay(x0,z0,param);
figure('Name','MUST - Matlab UltraSound Toolbox - Damien Garcia',...
    'NumberTitle','off');
subplot(334)
plot(dels*1e6,'.')
axis tight square
title('TX time delays in \mus')
xlabel('Element #')
subplot(3,3,[2 3 5 6 8 9])
plot((param.Nelements-1)*param.pitch/2*[-1 1],[0 0],'k','Linewidth',10)
axis equal ij
axis([-4e-2 4e-2 0 10e-2])
hold on
plot(x0,z0,'k*')
hold off
% Grid
x = linspace(-4e-2,4e-2,200);
z = linspace(param.pitch,10e-2,200);
[x,z] = meshgrid(x,z);
% RMS pressure field
P = pfield(x,z,dels,param);
% Final figure
imagesc(x(1,:)*1e2,z(:,1)*1e2,20*log10(P/max(abs(P(:)))))
caxis([-20 0])
c = colorbar;
c.YTickLabel{end} = '0 dB';
title('RMS pressure field')
colormap hot
hold on
plot((param.Nelements-1)*param.pitch/2*[-1 1]*1e2,[0 0],...
    'g','Linewidth',6)
axis equal ij off
axis([-4 4 0 10])
plot(x0*1e2,z0*1e2,'k*')
xlabel('[cm]')
hold off

end

